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From the molecular dynamics simulation of a system of hard-core disks in which an equilib- 
rium cell is connected with a nonequilibrium cell, it is confirmed that the pressure difference 
between two cells depends on the direction of the heat flux. From the boundary layer analysis, 
the velocity distribution function in the boundary layer is obtained. The agreement between 
the theoretical result and the numerical result is fairly good. 
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Although there has been a long history of nonequilib- 
rium statistical mechanics since Boltzmann introduced 
the Boltzmann equation, the understanding of nonequi- 
librium statistical mechanics is still in the primitive 
stages. ^^■^ The significant role of nonequilibrium physics 
in the mesoscopic region has been recently recognized. 
For example, there has been some important progress 
such as the Fluctuation Theorem^' ^ and the Jarzynski 
equality^ in mesoscopic nonequilibrium statistical me- 
chanics. In typical situations of mesoscopic physics, ma- 
terials are confined to narrow regions. Thus the boundary 
effects at the mesoscopic scale are important not only for 
nonequilibrium statistical mechanics but also for consid- 
eration of friction and lubrication.'' 

On the other hand, some macroscopic phenomenolo- 
gies for nonequilibrium steady states have been proposed. 
These are the Extended Thermodynamics (ET),^ the Ex- 
tended Irreversible Thermodynamics (EIT)^^^-'^ which is 
the combination of ET and information theory, and the 
Steady State Thermodynamics (SST).^^ It is interesting 
that both EIT and SST treat a common process in which 
a nonequilibrium cell is connected with an equilibrium 
cell"- 12 

The Knudsen effect is the phenomenon in which two 
equilibrium cells with different temperatures are con- 
nected by a small hole.^ The balance equation of the 
Knudsen effect is given by 

Pi P2 




(1) 



where Ti and Pi represent the temperature and the pres- 
sure in the cell, i, respectively. Although eq. (1) contains 
only macroscopic variables, it is in contrast to the ordi- 
nary thermodynamic balance condition where the pres- 
sures the two cells are equal. Such an exceptional condi- 
tion means that the mass balance is determined by the 
transportation of the gases in the small hole. Therefore, 
the relevance of predictions by macroscopic theory^"' 
for the generalized Knudsen effect in which an equilib- 
rium cell is connected with a nonequilibrium cell by a 
small hole is questionable. 

On the other hand, the explicit perturbative solution 
of the Boltzmann equation for hard spheres has been 
derived at the Burnett order of the heat flux.^^'^^ The 



Fig. 1. The schematic picture of the system of our simulation, 
where the right cell is at equilibrium with temperature T2 and 
the left cell is in a nonequilibrium state. The inset represents the 
detailed structure of the hole that has the gap 1.11. The unit 
of length in figures is the diameter of the rigid disk. We adopt 
periodic boundary condition in the y direction. 



quantitative accuracy and numerical stability of their so- 
lution in the bulk region has been confirmed by molec- 
ular dynamics (MD) simulation for hard-sphcres"'^^ and 
the extension to the tenth order shows that their sec- 
ond order solution is accurate even when the heat cur- 
rent is large. 1^ They have confirmed that the solution 
in the bulk region derived by information theory is not 
consistent with that of the Boltzmann equation. Kim 
and Hayakawa^^ have also discussed the nonequilibrium 
Knudsen effect. Their result denies the prediction of SST, 
but both theories predict that the osmosis AP defined 
by the AP = P"^"^ — P'^' with the xx component of the 
pressure tensor P"J^' in the nonequilibrium cell and the 
pressure P'^'' in the equilibrium cell is always positive 
regardless of the direction of heat flux. 

However, their simplification^^ using the bulk solution 
of the Boltzmann equation in the boundary layer is not 
acceptable. In fact, numerical simulations of the Boltz- 
mann equation of hard spheres^^ exhibit the discontinu- 
ity of velocity distribution function (VDF) in the bound- 
ary layer, but there is no discontinuity of VDF derived 
by Kim and Hayakawa. i'^' Furthermore, it is also well 
known that the gas temperature near the wall is differ- 
ent from the wall temperature,^''^' but both treat- 
ments^^'"'^^ ignore this fact. 

To clarify the truth of the nonequilibrium Knudsen 
effect, we employ the event-driven molecular dynamics 
simulation of hard-disks developed in refs. [20-23]. Let 
X and y be the Cartesian coordinates of the horizontal 
and the vertical directions, respectively. First we check 
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the validity of eq. (1) by MD. We adopt the diffusive 
reflection for the vertical walls away from the hole, the 
simple reflection rule for the wall between two cells and 
the periodic boundary condition for the horizontal wall. 
We connect two cells by a small hole, as illustrated in 
Fig. 1. For T1/T2 = 1.96, we have found the steady value 
of P1/F2 VW^i" = 0.9982 ± 0.0195, where the average 
area fraction and the number of the particle are 0.015 and 
10, 000 respectively. We determine the pressure based on 
the Virial theorem^'* and average in 10* collisions per 
particle. The stationary state is realized after 3 x 10^ 
collisions per particle and we use the data after that. 
We simulate the system until 10^ collisions have been 
performed per particle. 

Next, we simulate the nonequilibrium Knudsen effect. 
The number of hard disks and the average area fraction 
are the same values used to simulate the conventional 
Knudsen effect. We adopt the diffusive boundary condi- 
tion for the wall between two cells (see Fig. 1). We exam- 
ine the values of T1/T2 as 1.96 and 1/1.96. We cause both 
cells to divide into 10 equal parts, and we have confirmed 
that pressure based on the Virial theorem^"' is identical 
in each cell. This is consistent with the stationary con- 
dition.^^ In nonequilibrium cases the stationary state is 
realized after 10^ collisions per particle and we use the 
data after that. We simulate the system until 4.4 x 10^ 
collisions have been performed per particle. 

Our results of MD plotted in Fig. 2 indicate that the 
sign of AP(= — P'^'^) depends on the direction of 
energy flux. Actually, the stationary value AP in our 
simulation is given by AP/ < naTa >= 0.0279 ± 0.0184 
for T1/T2 = 1.96, and AP/ < nsTa >= -0.0152 ±0.0175 
for T1/T2 = 1/1.96, where < naTa > is the time average 
of naPa (Fig- 2). This result contrasts with the positive 
AP predicted by both SST^^ and Kim and Hayakawa. ^'^ 

Now let us compare the VDF of MD with the pertur- 
bative solution of the Boltzmann equation at the Burnett 
order obtained by Kim^'* for 2D hard disks. From now on, 
we restrict our interest to the data for heat flux J < 0, 
with T1/T2 = 1/1.96. As shown in Fig. 3, VDF obtained 
from MD in the bulk region of a nonequilibrium cell is 
almost identical to the theoretical VDF, where the VDF 
of MD is obtained from the average of particles exist- 
ing in ±1.05d from the center of nonequilibrium cell in 
the horizontal direction with hard disk diameter d. On 
the other hand, the VDF in the boundary layer in which 
we average the data of particles existing between 0.7c? 
and 2.8d apart from the right wall of the nonequilibrium 
cell deviates from the theoretical VDF,^* particularly for 
v,<0 (Fig. 4). 

Let us derive VDF in the boundary layer in the 
nonequilibrium cell. We assume that the VDF of the in- 
cident particles {v^ > 0) in the boundary layer is the 
same as the bulk distribution function {fNE)- On the 
other hand, we assume that the VDF for the particles 
reflected by the wall {v^ < 0) obey the Maxwellian. 

Thus, the distribution function Jbl in the boundary 
layer is given by 

r I X j fMB{n^o,T^o,^^) Vx<0 
\jNE\nx,Tx,J,^) Vx>Q, 
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Fig. 2. The time dependence of the pressure for (a) T1/T2 = 1.96 
and (b) T1/T2 = 1/1.96. Here P^x'^ in the nonequihbrium cell 
is plotted by plus points and its time average is shown by the 
dotted line. In each case the pressure P"'^ in the equilibrium cell 
is plotted by open square points and its time average is plotted 
by the dashed line. The theoretical P^x'' iu eq. (11) is shown by 
the solid line. For the normalization, we divide these values by 
< 712^2 >• 

(h-IO") 
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Fig. 3. Comparison of VDFs derived by simulation and Ref. [14] 
at the middle of the nonequilibrium cell. The plotted data are 
the substraction of Maxwellian from VDFs. 
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Fig. 4. Comparison of VDFs derived by simulation and Ref. [14] 
and Jbl in eq. (2) at the near the wall. The plotted data are the 
substraction of Maxwellian from VDFs. 
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with the mass of a hard-disk to, and 
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(4) 



with 6i = 1.03, 62 = 5.738 x 10"^ 63 = 4.946 x lO^^ 64 = 
4.313 X 10-4,65 = 3.452 x 10'^, be = 2.241 x 10"^ and 
Laguerre's bi-polynomial L1^{x).^'^ Since the heat flux is 
sufficiently smafl, we adopt the first order nonequihbrium 
VDF in the heat-fiux for Jne- 

There are three unknown variables, n^, n^, and in 
eqs. (3) and (4), while there are three relations, 



J = 



f 

/ d-VVx — fKLi-v)-, 



J dvvxfxLi^) = 0. 



(5) 
(6) 
(7) 



Here, the first two equations are definitions of the density 
n and the heat flux J, and the last equation represents 
the mass balance condition. Therefore, we can determine 
Hw, nx, and Tx from eqs. (5)-(7). The expansions in terms 
of J of the three variables become 

TO^/^ 

n + a ^,^ J, (8) 



J- 7 
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J, 



(9) 
(10) 



where a, 6 and c are constants to be determined. From 
eqs. (5), (6) and (7), we obtain the solutions of the linear 
simultaneous equations as (a, 6, c) ~ (0.32, —0.099, 0.84). 

The distribution function Jbl near the wall has 
thereby been determined by n, T^,, and J. Figure 4 is the 
comparison of Jbl with the results of MD. From Fig. 4, 
both fsL and the VDF from MD have a discontinuity 
at tij; = 0, as in the case of conventional boundary layer 
analysis.^* For reflective VDF {vx < 0), there is the small 
difference between the result of MD and the Maxwellian. 
We may deduce that it arises from collisions between 
particles because we measure the VDF a short distance 
from the wall. 

With the aid of fsL, AP can be calculated as 



APi 



0.412 J 



/7¥ 



= 0.415- 



J 



(11) 



From this equation and the value of the heat-flux in MD, 
we evaluate AP/ < nsTa >~ 0.0298 for T1/T2 = 1.96 
and AP/ < naTa >~ -0.0156 for T1/T2 = 1/1.96. The 
result agrees well that of the simulation. We can also 
obtain the temperature jump coefficient 7 ~ 0.97 which 
is defined through Tx=o — — jn^ m/Tyj. 

Now, let us discuss our result. There are some ad- 
vantages to employing MD as the numerical simulation. 
First, we can easily change boundary conditions for walls 
depending on our interest. The system of connecting two 
cells by a small hole is easily simulated by MD. Second, 
MD is suitable for high density simulation of gases. The 
nonequihbrium Knudsen effect for high density gases is 
an interesting subject for future discussions. We are also 
interested in the size effect of the hole on the transition 
from Knudsen balance to the thermodynamic balance. 
On the other hand, there are some disadvantages to em- 
ploying MD. Because of the small system size of our MD 
simulation, the fluctuation of the pressure is large. 

The advantage of our boundary layer analysis is that 
we can write the explicit form of VDF in the boundary 
layer in terms of the density, the heat-flux and the tem- 
perature of the wall. On the other hand, VDF by Sone et 
al^^ has an implicit form that is obtained as a numerical 
solution of an integral equation. The explicit VDF can 
be obtained by our simplification, but the validity of this 
method has not been confirmed. In fact, for 3D hard- 
sphere gases, our method predicts the temperature jump 
coefficient 7 ~ 0.72, while Sone et al^^ predicts 7 ~ 1.00. 
To check the validity of our boundary layer analysis, it 
will be necessary to employ 3D MD simulation for 3D. 

We also stress the difficulty of describing the nonequi- 
hbrium Knudsen effect by macroscopic phenomenology. 
In fact, our boundary layer analysis strongly depends on 
the boundary condition. 

In conclusion, the sign of AP depends on the direction 
of the heat flux. The approximated VDF in the boundary 
layer has been obtained from the assumption that the 
incident particles obey the bulk VDF and the reflected 
particles obey the Maxwellian. This agrees well with the 
simulation result. 
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